Vibrational properties of phonons in random binary alloys: An augmented space 

recursive technique in the k-representation 
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We present here an augmented space recursive technique in the k-representation which include 
diagonal, off-diagonal and the environmental disorder explicitly : an analytic , translationally invari- 
ant , multiple scattering theory for phonons in random binary alloys. We propose the augmented 
space recursion (ASR) as a computationally fast and accurate technique which will incorporate 
configuration fluctuations over a large local environment. We apply the formalism to Ni^Pd^s , 
Ni$$Cri2 and Ni^Pt^o alloys which is not a random choice. Numerical results on spectral func- 
tions, coherent structure factors, dispersion curves and disordered induced FWHM's are presented. 
Finally the results are compared with the recent itinerant coherent potential approximation (ICPA) 
and also with experiments. 
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I. INTRODUCTION 
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Over the years, there have been many attempts to 
develop adequate approximations for properties of ele- 
mentary excitations in disordered systems. Among these 
phonons are not only conceptually the simplest, but are 
also the most readily accessible to experimental verifi- 
cation. Neutron scattering experiments |l|, [2, y( have 
provided detailed information about lattice vibrations in 
random alloys. A satisfactory reliable theory is still lack- 
ing. The main reason for this is that, unlike the case 
of electrons in substitutionally disordered alloys, where 
the Hamiltonian can be expressed in a form where the 
disorder is diagonal in a real space representation (This 
is certainly true, for example, in a LMTO formalism 
in the absence of local lattice distortions), the disor- 
der in the dynamical matrix is essentially off-diagonal. 
To make things more complicated, the diagonal and off- 
diagonal disorders in the dynamical matrix are coupled 
by the force constant sum rule Qrr, — — J^r'^r^RR' > 
which ensures that no vibration can be excited in a uni- 
form translation of the crystal as a whole. This sum 
rule imposes an environmental disorder on the force con- 
stants. That is, the disorder in the diagonal element of 
the dynamical matrix depends upon its near neighbours 
or its immediate environment. Hence a reliable theory 
for phonon excitations will be that which includes all the 
three kinds of disorders explicitly. 

Let us look at the most successful mean field approxi- 
mation : the single-site coherent potential approximation 
(CPA), introduced first by Taylor 0. As the name it- 
self suggests, it is a single site approximation and per se 
cannot deal adequately with off-diagonal disorder. Sev- 
eral authors have proposed schemes for generalizing the 



CPA and their approaches include geometrically scaled 
off-diagonal disorder 0, linearly scaled off-diagonal 
disorder 0] and independent diagonal and off-diagonal 
disorder 0, HI EJ IHJ • Most of these schemes in practice 
lead to a single site CPA including off-diagonal disorder. 
These approaches suffer from two different kinds of draw- 
backs : first, there is no reason why the off-diagonal part 
of the dynamical matrix should scale either as the geo- 
metric mean or the arithmetic mean of the constituents. 
Secondly, often the extra assumptions lead to approx- 
imate Green functions which violate the essential Her- 
glotz analytic properties required to produce physically 
acceptable results. 

A complex function f(z) of a complex variable z is 
said to possess Herglotz analytic properties provided it 
satisfies the following properties 
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(i) The function is analytic everywhere in the complex z- 

plane, expect on the real z-axis. All its singularities 
therefore lie on the real z-axis. 

(ii) The sign of the imaginary part of f(z) is negative 

everywhere on the upper half of the z-plane and 
positive everywhere on the lower half of the z-plane. 

(iii) If the set of singularities of f(z) is bounded, then 

lim $tef(E) = G, E is real. 

z — E— >±oo 

The augmented space approach suggested by Mook- 
erjee |l2j provided a very interesting starting point for 
the generation of appropriate approximations. The 2- 
site CPA proposed by Yussouff and Mookerjee for 
model systems and subsequently generalized by Mooker- 
jee and Singh 0,0] to realistic alloys was one success- 
ful approach. While it retained the Herglotz properties 
of the approximate Green functions, the generalization of 
the 2CPA to larger clusters violated the lattice transla- 
tional symmetry of the configurationally averaged Green 
function for homogeneous disorder. This was overcome 
subsequently by the traveling cluster approximation of 
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Kaplan and Gray |8| also based on the augmented space 
method. Recently, Ghosh et.al have proposed a near- 
est neighbour traveling cluster CPA and have applied it 
to phonons in NiPt and NiPd alloys. In this communi- 
cation we shall propose a different approximation proce- 
dure. We shall start from the augmented space method 
and use the recursion method of Haydock et.al [2(| to 
obtain the configurationally averaged Green functions. 
The termination of the continued fraction expansion will 
constitute the approximation. This will not only retain 
the Herglotz analytic properties of the approximate aver- 
aged Green function, but also include the effect at a site 
of its neighbourhood, the size of which we can control. 
We shall incorporate the effect of the very distant envi- 
ronment by the use of accurate termination schemes pro- 
posed e.g. by Haydock [2l| , Luchini and Nex |22( or Beer 
and Pettifor [2j|. Since we shall incorporate the lattice 
translation symmetry in the full augmented space (which 
is characteristic of homogeneous disorder) |24j within our 
approach, the drawback of the original cluster-CPAs used 
by Singh and Mookerjee [lH. Il4| will be overcome. Fur- 
ther, we shall use the local point group symmetries of 
the lattice and the configurations on it to drastically re- 
duce the rank of the Hilbert space on which the recur- 
sion takes place (see Saha et.al ^3)- This will allow 
us to accurately account for large environments around 
a particular site. One of the strengths of the proposed 
method which will represent a major step forward in the 
theory is the possibility of including random fluctuations 
in force constants beyond the nearest neighbours. While 
in certain representations the Hamiltonian of electronic 
systems can be seen to be short-ranged, this is not so for 
dynamical matrices. The recursion method in augmented 
space can include beyond nearest neighbour randomness 
in force constants without much computational expense. 
In our work on NiPt and NiPd we have extended disor- 
der upto second nearest neighbours to illustrate this. It 
is not immediately clear how easy it would be to extend 
the method proposed by Ghosh et. al to larger sized clus- 
ters. We propose the ASR as a computationally fast and 
accurate technique which will incorporate configuration 
fluctuations over a large local environment. 
In Section II, we shall introduce the basic formalism. 
In Sections III we shall present results for Ni^Pd^, 
Nig,$Cri2 and Ni^Pt^ alloys and compare them with 
experiment. The choice of the systems is deliberate : 
NiPd has predominantly mass disorder, NiCr predomi- 
nant disorder in the dynamical matrix ; NiPt has large 
disorder both in the mass and the dynamical matrix. 
Concluding remarks are presented in Section IV. 



II. FORMALISM 

A. The augmented space formalism for phonons 

The basic problem in the theory of phonons is to solve 
a secular equation of the form : 

(Mw 2 -D) u(R,lj) = 

where u a {R,uj) is the Fourier transform of u a (R,t), the 
displacement of an atom from its equilibrium position R 
on the lattice, in the direction a at time t. M is the mass 
operator, diagonal in real space and D is the dynamical 
matrix operator whose tight-binding representation is of 
the form : 

M = J2 m R <W P R C 1 ) 

R 

D = E *& Pr + E E *& T ™ ( 2 ) 

R R R'^R 

along with the sum rule: 

<R = E <R> ( 3 ) 

R'^R 

Here Pr is the projection operator \R)(R\ and Tr_r< is 
the transfer operator \R)(R'\ in the Hilbert space H 
spanned by the tight-binding basis R,R' specify 

the lattice sites and a, (3 the Cartesian directions. m,R is 
the mass of an atom occupying the position R and 
is the force constant tensor. 

We shall be interested in calculating the displacement- 
displacement Green function in the frequency-wavevector 
space , which in the absence of disorder in the system has 
the diagonal element 

G(k, k', oj 2 ) = G(k, oj 2 )6(k - k') 
and for the present case 

G(k,w 2 ) = (k| (Mw 2 - D) _1 |k) 
where |k) is a state in the reciprocal space given by : 

Ik) = -4=yexp(-tk.R)|i2) 

Since the mass matrix M is perfectly diagonal, we can 
write 

G(kV) = (klM-Wo^I-M-^DM- 1 / 2 )" 1 

...M~ 1/2 |k) (4) 

where 

M- 1 /a=^ m -V2^ Pfl 

R 
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The equation J2J looks exactly like the Green function for 
the electronic case with M _1 / 2 DM -1 / 2 playing the role 
of Hamiltonian H, uj 2 in place of energy and M _1 / 2 |k) is 
the starting state of recursion. 

Let us now consider a binary alloy A x B y consisting 
of two kinds of atoms A and B of masses ttia and 
tub randomly occupying each lattice sites. We wish 
to calculate the configuration-averaged Green function 
<§; G(k,uj 2 ) 3>. We shall use the augmented space for- 
malism (ASF) to do so. Since the disorder is homoge- 
neous, averaged <C G(k, uj 2 ) 3> is also diagonal in re- 
ciprocal space representation 123. The ASF has been 
described in great detail earlier [23. We shall indicate 
the main operational results and refer the reader to the 
above monograph for further details. The first operation 
is to represent the random parts of the secular equation 
in terms of a random set of local variables {n R } which 
are 1 if the site R is occupied by an A atom and if 
it is occupied by B. The probability densities of these 
variables may be written as : 

Pr(n R ) = x 8(n R - 1) + y S(n R ) 

= (-1/tt) SMThI (nnl-Nn)- 1 \U) (5) 

where x and y are the concentrations of the constituents 
A and B with x + y = 1. N R is an operator defined on 
the configuration space <p R of the variable n R . This is of 
rank 2 and is spanned by the states {|tfl), \ Ir)} '■ 



x p r + yp l R 



xy(r R l 



Let us now carry out the ASF operations in some de- 
tail. The mass m R for site R can then be expressed as 

-1/2 -1/2 — 1/2 X 

m R = m A n R + m B (l-n R ) 



A 
,-1/2 



(8m) 



-1/2 



llR 



where 



Therefore 



(8m) 



-1/2 



m 



-1/2 



m 



-1/2 



R 



\m B 1 ^ 2 + n R (8m) x l 2 



8 a p Pr 



(6) 



In augmented space formalism, in order to obtain the 
configuration average we simply replace the random vari- 
ables n R by the corresponding operator N R associated 
with its probability density (as in equation |SJ) and take 
the matrix element of the resulting operator between 
the reference states. The justification is sketched in the 
Appendix A. For a full mathematical proof the reader is 
referred to [12 . 



Using the above in equation JBJ) we get, 



MT 1/2 = m~ 1/2 I®I- 



,-V2 V- i 



Y.PR® p R 



in 



-1/2 



E Tr®Pr 



(7) 



R 



where 



-1/2 -1/2 
m x = x m A + y m B 

m 2 = (y — x) ((5m) -1 / 2 

77i 3 1 ^ 2 = ^fxy ((5m) -1 / 2 
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Similarly the random off-diagonal force constants ^ a RR , 
between the sites R and R' can be written as : 



$ a/3 
^RR' 



^f A n R n w + $f B (l - n R )(l - n R ,) 
... + [ n R (l - n Rr ) + n R , (1 - n R ) 



■■■ + (<b-<b) (n R + n w ) 



a/3 



S«/3 



(8) 



Let us define the following 



a/3 

(i) 

a/3 
(2) 



x*fA-yn P B + (y-xm p B 



a/3 



a/3 
AA 



af3 
BB 



2$ 



a/3 
AB 



In augmented space the off-diagonal force constant ma- 
trix becomes an operator : 



a/3 _ 



E 

RR' 



+ ^{(y- X )(p R +P R ,) 
+ y/xy(T^+T^)}+^ ) {(y-x) 2 p R p R ,+ 
xy(y - x) [p R 7# +p R , 7^) + xy t£t#} 



T RR > 



E K%^T RR , 



RR' 



The sum rule, 



(9) 



a/3 
RR 



Z < ^RR' 
R'^R 



gives the diagonal element of the dynamical matrix 



n R 



■ N R =x p R 
= x I 



■yp l R 



^y(r R l 



(y -x)p l R + ^y~ T R 



T R > 

n 




E *& 



(10) 
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The total dynamical matrix in the augmented space is : 

6 = -E I E n%\®p«. + E C®^ 

(11) 



The augmented space theorem |l2j now states that the 
configuration averaged Green function <C G(k, w 2 ) 3> 
may be written as : 



«G(kV)» = J G(k,^ 2 ,{n fl }) Yl Pr(n R )dn R 
= (k®{0}| G(kV,{^})|k®{0}) 

= <k ® {0}| NT 1/2 (w 2 I - M- 1/2 DM- 1/2 ) ~* 1VT 1/2 |k ® {0}) 

I 



(12) 



where M -1 / 2 and D are the operators which are con- 
structed out of M -1 / 2 and D by replacing all the random 
variables n R (or n^' ) by the corresponding operators N R 
(or Nri) as given by equation Q and These are the 
operators in the augmented space = W ® <&. The state 
|k(E>{0}) is actually an augmented space state which is the 
direct product of the Hilbert space basis and the Config- 
uration space basis. The configuration space <fr = Yl% 4>R 
is of rank 2 N for a system of N-lattice sites with binary 
distribution. A basis in this space is denoted by the car- 
dinality sequence {C} = {Ri, R2, ■ ■ ■ , R c } which gives us 
the positions where we have a | J.) configuration. The 
configuration {0} refers to a null cardinality sequence i.e. 
one in which we have | f) at all sites. 

Using the operator representation for M" 1 / 2 we get : 



M- 1 / 2 |k®{0}) 



,-1/2 



{0}) 



,-l/ 2 



{R}) = \i} 



where a configuration state is denoted by its cardinality 
sequence {G}. We have also used the short hand notation 



^L^TexpHk-i?) \R ®{G}) 



{C}) 



The ket |1} is not normalized. A normalized ket |1) is 
given by 



ID = 



I 1 } 

VI 

With the definitions 



im=(-)" /2 " «>+(^r /2 11™ 

1 Ml \ m / \ m J 



we may rewrite equation l|12|) as : 



«G(k,c 2 )» = (1| (w ,2 I--D eff )- 1 |1) 



(13) 



where, us' 2 = moj 2 and D e // = I D X 

This equation is now exactly in the form in which 
recursion method may be applied. At this point 
we note that the above expression for the averaged 
•C G(k, u}' 2 ) ^> is exact. The recursion transforms the 
basis through a three term recurrence relation as : 



\<h) = |i) ; I0o> = 

|^n+l) = D e //|0n) - - 6*|0 n -l) (14) 

The averaged Green's function can then be written as a 
continued fraction : 



<G(k,w 2 ) >= 



63 



/2 



/2 

W — (22 — 



- ajv -r(k, w ' 2 ) 



(15) 

where r(k, uj 12 ) is the asymptotic part of the continued 
fraction, and 



< (1/m) > 1 = m 



1 = 



I®I- 



-1/2 



E 



a„(k) = 

&n0<) = 



((/) n |D e// |0„ 
(0n|0n) 

(0n-l|0n-l) 



and 



6? = 1 



(16) 
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To implement the above recursion, we require to know 
the effect of the operator D R f / on a general state in aug- 
mented reciprocal space |3l| . Some of the main opera- 
tions are shown in Appendix B. 

So far the expression for the averaged Green function 
is exact. Approximations are introduced at this stage 
for its actual numerical evaluation. The mean-field the- 
ories essentially obtain the self-energies because of dis- 
order scattering, self-consistently and approximately and 
then calculate the averaged Green function either from 
the Green function without disorder or the virtual crys- 
tal Green function. The Coherent potential approxima- 
tion (CPA) proposed by Soven 1291 . the cluster CPA 
proposed by Mookerjee and Singh [l^L IT3 |. the travel- 
ing cluster approximation (TCA) proposed by Mills and 
Ratanavaraksa 9] and Kaplan et. al [10| and the itinerant 
CPA (ICPA) proposed by Ghosh et.al [l7| basically all 
belong to this category. The latest work referenced rep- 
resent the most sophisticated version of the mean-field 
theories. We shall propose an approximation that will 
start from the infinite continued fraction and approxi- 
mate its asymptotic part by an analytic termination pro- 
cedure. The coefficients a n , b n are calculated exactly up 
to a finite number of steps and the asymptotic part is 
then replaced by a terminator : r(k, uj' 2 ) ~ T(k, cj' 2 ). 
The concept of terminators is described in the appendix 
C, where further details of estimating T(k, oj' 2 ) are also 
provided. Haydock and co-workers |2l| have carried out 
extensive studies of the errors involved and precise esti- 
mates are available in the literature. Several terminators 
are available and we have chosen to use that of Luchini 
and Nex If we calculate the coefficients up to the 
n-th step exactly the first 2n moments of the density of 
states are reproduced exactly. The terminator ensures 
that the approximate Green function has Herglotz ana- 
lytic properties which also tells us that the approximate 
density of states is always positive definite and the spec- 
trum is always real. The terminator is also so chosen that 
the asymptotic moments are also accurately reproduced. 
This is a generalization of the method of moments, with 
the additional restriction that the asymptotically large 
moments are also accurately obtained. 
In the absence of disorder in the problem, the Green func- 
tion for a given mode is of the form : 



G (kV 2 )= ,., 



-c 2 (k) 

The spectral function A (k, a/ 2 ) is a delta function of 
the form 8{u>' 2 — Wg(k)). If we write : 



E(k,a/ 2 ) = ai (k)-co 2 (k) 



&l(k) 



lu - a2(k) 



= ai (k)-a; 2 (k) + ( T(k, W ' 2 ) 



(17) 



< G(k,u' 2 ) > = G (k,cy 2 - £(k,cy 2 )) 

Obviously from above £(k, a/ 2 ) is the disorder induced 
self-energy. Damped vibrations occur with reduced fre- 
quencies at w (k) which are the solutions of the implicit 
equation : 



w 2 (k) - ai(k) - 5Re o-(k,w 2 (k)) = 



and their disorder induced widths are : 



W(k, Wo 2 (k)) 



3m a(k,w 2 (k)) (18) 



The average spectral function <C A\ (k, w 2 ) ^S> for a 
mode labelled A is related to the averaged Green function 
in reciprocal space as : 



<A A (kV 2 )>=-i lim hm\<z:Gx(k,Lu' 2 -i5)-»\ 

7T ,5—0+ L I J 

(19) 

The averaged density of states is given by : 



< n{J) >= 



2uj' 

IT 



E 



— < A x (k,uj' ) > 



BZ 



Then, 



Here A labels the particular normal mode branch and BZ 
is the Brillouin zone. 

The dispersion curves for different modes are then ob- 
tained by numerically calculating the peak frequencies of 
the spectral function. This averaged spectral function 
gives, in principle, a proper description of the dynam- 
ics but it does not involve any weighting by scattering 
lengths. The dispersion curves so obtained are nearly 
the same as those obtained experimentally from the peak 
frequencies of the coherent structure factors S co h- This 
is because the coherent structure factors are nothing but 
the averaged Green functions weighted by the coherent 
scattering lengths. The intensities and the line shapes 
measured from S co h and the imaginary part of Green 
function may differ significantly, but the peak positions 
will generally differ little. 



B. The coherent scattering structure- factors 

Experimental determination of the phonon dispersion 
and line-widths are deduced from the averaged coherent 
scattering structure factors. The expression for these can 
be written as : 



< S coh (k,uj 2 ) >= 3m«bG(q, w 2 )b» (20) 

7T 



6 



here, thermal neutrons with wave- vector k gets scattered 
to final state of wave- vector k', q = k-k'+Q with Q 
being a reciprocal lattice vector. The energy lost by 
the incoming neutrons are taken up by the phonons : 
(h 2 /2M n )(k 2 - k' 2 ) = huj. f(w) is the Bosc distribu- 
tion function and 

b = ^2 b R 5 afj Pr 

R 

where 6_r is the scattering length of the nucleus occu- 



pying the site R. This is a random variable taking two 
values b A or b B depending on which kind of atom sits at 
the site labelled R. For comparison with experiment we 
have to calculate —(1/ir) 3m <C (b G(q,w) b) 
rather than the spectral function. For ordered materials 
the two are proportional, but if the scattering lengths 
are themselves random then although this has very little 
effect on the dispersion curves, it does affect the line- 
shapes and line widths. We can easily implement such 
an average within the ASF : 



<C b G(q,w 2 ) b > = (q< 



{0}| b [Mlo 2 
I 



-l _ 
b |q< 



{0}> 



where 

b = <b» / ® I + (y - x){b B - b A )J2P l R® P R 

R 

... + ^y(b B - b A ) r R ® p R 

R 

Carrying out algebra similar to the one for the averaged 
Green function, we obtain : 



where 



bG(q,c 2 ) b»= (l b | (u^I-D^) 1 |1 6 ) (21) 



( y) II{0})+ (t) 11 {R}) 



with, 



X x 1/2 = x m A 1/2 b A + y m B 1/2 b B 
X^ 1/2 = (y - x) (m^ 1/2 b A - m B 1/2 b B 



x 3 1/2 = V%y 



(m A 1/2 b A 



m 



-l/2i 



Also 



where, 



b 2 b 2 b 2 

A = < — >=x^-+y-§- 

m m A m B 



Xuj 2 and D b eff = W D W 



(22) 



W 



mi 
X 



-1/2 



-1/2 



m 2 
X 



-1/2 



Tl 
R 



(23) 



The subsequent recursion calculation follows the iden- 
tical steps as for the averaged spectral functions. We 
have chosen a second neighbour force constant model, 
with dynamical matrices fitted to reproduce the disper- 
sion curves. The disorder induced widths are the quanti- 
ties which are more sensitive to the effect of randomness 
as compared to the frequencies (i.e. dispersion curves ), 
and as such will be one of the focus of this work. In or- 
der to extract the full width at half maxima ( FWHM 
) , we have fitted the coherent structure factors to Loren- 
zians exactly as the experimentalists do to extract the 
same. The advantage of including the scattering length 
fluctuation will be clear when we will show the nature 
of the line widths for Ni^Pd^ alloy with and without 
inclusion of the scattering length fluctuation (i.e. Calcu- 
lating the widths once by fitting the spectral functions 
to Lorcnzians and then the structure factors to Loren- 
zians ). Our aim in this communication is to propose 
the augmented space recursion as a useful technique to 
study effects of diagonal, off-diagonal and environmental 
disorder. Accurate model building or obtaining the force 
constants from first-principles total energy calculations 
will be postponed for future work. 



In the following three Subsections, we present our cal- 
culations on Ni 55 Pd45, iVi 88 Cri2 and Ni 50 Pt 50 alloys. 
The choice is not arbitrary. Mass disorder dominates in 
NiPd alloys while force constant disorder is large in NiCr 
alloys. NiPt alloys have large disorder both in mass and 
force constants. Since in the phonon problem we have 
both kinds of disorder, it would be interesting to note 
the interplay between them in this series of alloys. The 
concentrations are chosen so that we may compare our 
results with existing work. 
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FIG. 1: Total coherent structure factors in different directions with different branches for NissPd45 . In each of the different 
directions and branches, the various curves indicate the total structure factors for various £ values starting from the lowest 
value to the edge of the Brillouin zone. In [£ 0] direction Tl and T2 modes are degenerate, in [f £ 0] direction L and Tl 
modes are degenerate and in [£ ( £] direction all the three modes are degenerate. The y-axis is in an arbitrary scale with 
heights scaled to the maximum height. Different curves for different £ values are shifted along the x-axis in order to facilitate 
vision. 





Ni 


Pd 


Pt 


Cr 


Atomic number 


28 


46 


78 


24 


Atomic mass (amu) 


58.71 


106.4 


195.09 


51.996 


Free atom valence config. 


3d 8 4s 2 


4d w 


5d"6s L 


3d a 4s L 


Lattice constant (fee) 


3.524 


3.8904 


3.924 


3.68(fcc) 
2.89(bcc) 


Elastic constants at 
296° K (I0 12 dyne /cm) 










Cn 


2.461 


2.270 


3.467 


3.5 


Cl2 


1.501 


1.759 


2.507 


0.678 


C44 


1.220 


0.717 


0.765 


1.010 


n-n force constants 
(in units of dyne/cm) 










IXX 


17319 


19337 


26358 


37483 


IXY 


19100 


22423 


30317 


17453 


IZZ 


-436 


-2832 


-7040 


-13229 


n-n-n force constants 
(in units of dyne/cm) 










2XX 


1044 


1424 


4926 




2YY 


-780 


210 


-537 





TABLE I: General properties of fee Ni, Pd, Pt and Cr. The 
force constants for Ni, Pd and Pt are taken from |27| and that 
for Cr is taken from 14] 



III. RESULTS AND DISCUSSION 

A. The Ni 55 Pd 4 5 Alloy 

A look at Table [I] immediately shows us that for NiPd 
alloys, the dominant disorder is in the mass. Force con- 
stants in Pd are only about 15% larger than those in Ni. 
We shall choose the Ni5sPd45 alloy for the application of 
our formalism developed in section II. This particular al- 
loy has already been studied within the ICPA by Ghosh 
et. al , CCPA by Mookerjee and Singh Q and exper- 
imentally by inelastic neutron scattering. The alloy 
forms a continuous series of face centered cubic solid so- 
lutions of all concentrations and there are no indications 
of long range order down to 0°C. 

In Figure ^ we display the coherent scattering struc- 
ture factors obtained from our recursion calculation along 
the highest symmetry directions ( [C00], [CC0], [CCC] )> 
£ = |fc|/|fc maa; | for different branches. For a particular 
direction and branch the different curves indicates the 
spectral functions for various C points starting from the 
lowest value ( i.e. £ = ) to the edge of the Brillouin 
zone ( i.e. £ = 1 in units of 2n/a ). The first thing 
to note is that the structure factors are (in contrast to 
Lorentzian shape) often asymmetric near the resonances. 
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FIG. 2: Dispersion Curves ( frequency v vs. reduced wave 
vector £ ) for Ni5sPd45 calculated from recursion (solid line). 
The force constants used are given in the text. The filled 
circles are the experimental data Q . In all the three panels 
the thin dotted lines span the FWHM's . 



The asymmetries can be described as a tendency of more 
scattering to occur near the resonance frequencies. In 
other words the shape of a mode with a frequency slightly 
lower or higher than that of a resonance tends to have 
a second peak or wide tail over the resonance region. 
In fact if one looks at the [CCO] L or Tl ( doubly de- 
generate ) and [CCC] L or Tl or T2 (3-fold degenerate) 
branches, the shape of a doubly peaked structure factor 
is much more clear. Out of these two peaks, one peak 
corresponds to the dispersion curve for the longitudinal 
mode (L) and the other peak to the transverse mode (T). 
Experimentally, for some neutron groups corresponding 
to transverse phonons with frequencies just below the 
lower resonance, definite asymmetry to the right was ob- 
served. Such asymmetries are clearly observed for the [£ 
00]T and [££0]T1 branches. It is important to note that 
the structure factors have a pronounced k and branch 
dependence. 

In figure|21we display the dispersion curves, which were 
constructed by numerically determining the peaks in the 
coherent scattering structure factor. In this communica- 
tion our main focus is the development of the augmented 
space recursion method. Accurate determination of the 
force constants shall be left for the future. Ghosh et.al 
|17| have attempted much more detailed determination 
of the force constants. For the time being we have used 
the same parametrization of the force constants as they 
did. These dispersion curves ( solid lines ) are compared 
with the experimental results Q ( filled circles ). The 
dotted lines span the calculated FWHM's. The proce- 
dure of calculating FWHM's has already been discussed. 
The asymmetry in the widths is again clearly observed in 
the two transverse branches quoted above. The results 
are in good agreement with the experiment for all the 
three symmetry directions and for each branch. Agree- 
ment can be achieved by varying only one parameter in 



the force constant matrix. This suggests that the force 
constant disorder is weak and the system is dominated 
by the mass disorder, as is clear from the numerical val- 
ues of the parameters given in Table [I] . If one looks at 
the previous results for the dispersion curves ( i.e. VCA, 
CPA and ICPA curves ) |17|. it will be clear that in the 
low wave vector regime, there is no distinction between 
these results and ours, Because the self averaging of both 
mass and force constants over a single wavelength re- 
duces the CPA, ICPA and the ASR results to become 
close to the VCA. However as we move toward high wave 
vectors, the VCA curve deviates from the experimentally 
observed one and lies lower in frequency as compared to 
the ICPA and the ASR results. The reason is that VCA 
uses an averaged mass. In contrast to this, for those the- 
ories which capture the effect of mass fluctuation (as do 
the ICPA and the ASR), the lighter atoms (Ni in this 
case) dominate in the high wave vector region and push 
the frequencies up. That is why our results agree very 
well across the Brillouin zone. 

The FWHM's are much more sensitive to approxima- 
tions as compared to the frequencies. These are shown 
in figure ©. The FWHM's shown in the left arc those 
which have been calculated without including any scat- 
tering length fluctuation, while in the right are those 
where the fluctuation has been included. The circles 
along with the error bars are the experimental data pj. 
It is obvious that the nature of the line widths are not the 
same in the two cases, rather the one including scattering 
length fluctuation is matching more closely with the ex- 
perimental data than the one without including the fluc- 
tuation. That should be obvious because the experimen- 
talists do include this fluctuation. Our results show very 
strong branch and wave-vector dependent widths and in 
good agreement with the experimental results of Kami- 
takahara and Brockhouse except in the [CCC]L mode. The 
reason for this may be because of the highly asymmet- 
ric line shapes in the [CCC]L mode. The single site CPA 
yields branch and k independent widths. It cannot cap- 
ture the essentially off-diagonal disorder of the problem. 
The ICPA and the ASR manages to capture this fea- 
ture. One should note that the structure factors are often 
asymmetric in shape and the usual Lorenzian fits carried 
out by most people may not be valid. 



B. The Nig 8 Cr 12 alloy 

We shall choose this alloy as being the nearest to that 
studied experimentally by Bosi et.al psj. Determination 
of the force constant matrices for this alloy becomes diffi- 
cult, because pure Cr is body centered cubic, but alloyed 
with Ni, up to 30% Cr it forms face centered cubic alloys. 
The force constants of pure Cr may be nothing like those 
of Cr in this alloy. Until we are able to determine these 
from a more first-principles type approach, our determi- 
nation of the force-constants for this alloy will remain 
tentative. We shall consider a hypothetical fee Cr, whose 
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FIG. 3: Full widths at half-maximum for the NiPd alloy as function of frequency for different directions in k-space and different 
modes without (left) and with (right) the inclusion of the scattering length fluctuation. The filled circles along with the error 
bars are the experimental data |IJ 



force constants are related to the elastic constants of bcc 
Cr via : 



C n + C12 = 4(/, - f t - - ft) /a 
C n -C 12 = (fi + 5f t , + f t )/a 
C44 = (fi + ft>+2ft)/a 

The values of C n, C12 and C44 are taken from 
Leibfried and Breuer [29j (given in Table I). 
It has been observed that the spectral functions and 
the structure factors for Ni 88 Cri 2 has strong evidence of 
branch dependent widths as also asymmetry in certain 
directions. This lends credence to our belief that force- 
constant disorder leads to both asymmetry and strong 
wave- vector and frequency dependence of the line-shapes. 

The influence of force constant disorder may be demon- 
strated more prominently in the dispersion curves and 
widths. In figure^ we display the dispersion curves along 
with the FWHM's using the force constants of Table [I] 
The procedure has already been discussed in the previ- 
ous section. These dispersion curves compare well with 
the experimental results |2a as well as the 2 CPA results 
of Mookerjee and Singh ([15|) (filled circles). The dotted 
lines span the calculate FWHM's. It should be noted 
that in the low frequency region, the widths are small 
but start to become significant as the phonon frequency 
increases. The widths are comparatively larger in the 
[COO] L, [CCC] L and [ CCO] T2 branches for high (-values. 
Looking at the dispersion curves, one should notice that 
the behaviour of the natural widths were somehow com- 
plemented in the behaviour of the frequencies. There is 
little evidence of resonances. This is expected, since clear 



cut resonances are characteristics of large mass disorders 
only. 

In figure [S] we show the FWHM as a function of fre- 
quency. It is clear that there is strong evidence of mode 
and k-dependencc. The FWHMs are very large and 
asymmetric for the longitudinal modes near the band 
edge frequencies. 

It is obvious from the above discussions that the force 
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FIG. 4: Dispersion Curves ( frequency v vs. reduced wave 
vector £ ) for NiggCri2 calculated from recursion (solid line). 
The force constants used ar e gi ven in the text. The filled 
circles are the 2 CPA results [lql . In all the three panels the 
thin dotted lines span the FWHM's . 
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FIG. 5: Full widths at half-maximum for the NiCr alloy as 
function of frequency for different directions in k-space and 
different modes. 



constant disorder plays a significant role in Ni 88 Cr 12 ; 
and a theory capturing only mass disorder effect ( e.g. 
like CPA ) fails to provide various essential features. 



C. The NisoPtso alloy 

Being encouraged by the right trend of theoretical re- 
sults toward the experimental results in the NissPd45 
and NiggCr^ alloys, where either of the two disorders 
- diagonal and off-diagonal dominates, we now apply our 
formulation to NiPt alloys where both disorders are pre- 
dominant. The mass ratio mp t /nr/vi is 3.3 (quite large 
compared to previous alloys) and the nearest neighbour 
force constants of Pt are on an average 55% larger than 
those in Ni. Tsunoda et.al Q have studied this sys- 
tem thoroughly covering a wide range of concentration 
(x=0.05, a?=0.3, x=0.5) by inelastic neutron scattering. 
In our case, we have considered x=0.5 where we expect 
the disorder induced scattering to have the strongest ef- 
fect. 

In this case, the spectral functions as well as the struc- 
ture factors show few extra features : Even in [(0 ]L, 
[( ]T and [( ( 0]T2 modes, unlike the previous two 
cases both the functions have one usual well defined peak 
( observed more clearly in the middle-regime of the Bril- 
louin zone ) along with a weakly defined peak with no 
gap in between. The occurrence of such a weakly defined 
peak is due to the inclusion of force constant disorder. 
Ghosh et.al have argued that it is entirely because 

of the off-diagonal disorder in the force constants. We re- 
fer the reader to their paper for the detailed arguments. 
Here we note that the feature is equally well reproduced 
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FIG. 6: Dispersion Curves ( frequency v vs. reduced wave 
vector ( ) for NisoPtso calculated from recursion . The force 
constants used are given in the text. The solid lines are the 
L-branch in all the three panels, the dashed lines are the T- 
branch in the left and right panels. In the [( ( 0] direction 
the dashed line indicate the Tl branch while the dot-dashed 
line indicate the T2 branch. The filled circles are the ICPA 
results 113] • In all the three panels the thin dotted lines span 
the FWHM's. 



in our augmented space recursive technique as well, why 
this should be so ? The effect of force constant disorder 
can be understood more clearly by looking at the disper- 
sion curves and widths. 

In figures El and d we display the dispersion curves and 
widths respectively obtained in the recursion using the 
force constants as given in the text. The solid lines are 
the L-branch in all the three panels, the dashed lines are 
the T-branch in the left and right panels, In the middle 
panel the dashed line indicate the Tl-branch while the 
dot-dashed line indicate theT2-branch. We have used 
the same parametrization of the force constants as used 
by Ghosh et al. These dispersion curves (solid lines) are 
compared with that calculated in the ICPA (filled circles) 
by Ghosh et. al ^3 • The aim of this work was to establish 
the ASR as a computationally fast and accurate method 
for phonon calculations for random alloys. Our ultimate 
goal is to obtain the force constants from first principles 
calculations. 

For all the three panels the thin dotted lines indicate 
the FWHM's. Unlike the previous two cases, the disper- 
sion curves in this case have very different characteris- 
tic features. The splitting of the curves in all the three 
symmetry directions is the main feature. This is due 
to strong resonances, a feature of large mass disorder. 
Also as one can see that near the resonances (around 
4 THz) the FWHM becomes very large. Tsunoda finds 
this resonance near 3.8 THz, while the CPA gives a rather 
lower value of 3 THz. In addition to these features , It 
has been observed that around 7 THz, the structure fac- 
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attempt to obtain the dynamical matrix itself from more 
microscopic theories. 
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APPENDIX A: THE AUGMENTED SPACE 
FORMALISM 

Let f{n R ) be a function of a random variable tlr , 
whose binary probability density is given by : 

Pr{n R ) = x A S(n R ) + x B S(n R - 1) 
We may then write : 



FIG. 7: Full widths at half-maximum for the NiPt alloy as 
function of frequency for different directions in k-space and 
different modes including the scattering length fluctuation. 



tor has a third small peak split from the main branch. 
This evidence of a weak resonance was also speculated 
by Mookerjee and Singh (H3)- The overall agreement of 
our dispersion curves with those calculated in the ICPA 
is good . 



IV. DISCUSSION AND CONCLUSIONS 

We have set up the augmented space recursion in re- 
ciprocal space for the study of phonon dispersion and 
disorder induced line-widths and line-shapes for random 
binary alloys. The technique takes into account both di- 
agonal disorder in the masses, the off-diagonal disorder 
in the force constants and the environmental disorder in 
the diagonal term of the force constants arising out of 
the sum rule. The approximation involving termination 
of continued fraction expansions of the Green function 
retains the essential Herglotz analytic properties. We 
have applied the method to three classes of alloys : NiPd 
where mass disorder dominates, NiCr where force con- 
stant disorder dominates and NiPt where both dominate. 
Wherever possible we have compared our results with 
neutron scattering data as well as the most sophisticated 
mean field theory recently proposed by Ghosh et.al (H3)- 
Both qualitatively and quantitatively our results agree 
well with the available data. We propose the technique 
as a computationally efficient method for the study of 
phonons in disordered systems. Our approach here made 
no attempt to obtain the force constant themselves from 
first principles, but rather resorted, as others did earlier, 
to fitting them from experimental data on the constituent 
metals. Our future endeavor would be to rectify this, and 



Pr(n R ) 



-— 3m(f_R \(n R I — Nr)~ x \ *\ r ) 

7T 



Here, the operator Nr acts on a space spanned by 
the eigenvectors |0) and |1) of Nr, corresponding to 
eigenvalues and 1 ; | 1r) = ^/xa\G) + ^/xb\1) is 
called the reference state. Its orthogonal counterpart is 
I Ir) = \/x~b\0) — y/xX\ 1) • The representation of Nr in 
this new basis is : 



Nn 



Now, 



< f(n R ) > 



f(n R )Pr(n R )dn R 



= — c$m f(n R )(U\(n R I-N R )-> 

* J-oo 

• ■ ■ I ] R )dn R 
= -isjfm E E /("*XtH |A> 

n A=0,1 A'=0,l J -°° 

...(A|(n H /-iV fl )- 1 |A')<A'| U)dn R 
= E (tfl |A)/(A)(A| U) 



A=o,r 

= (Tfl|f|Tfl) 



(Al) 



Here f is an operator built out of f{n R ) by simply 
replacing the variable n R by the associated operator N R . 
The above expression shows that the average is obtained 
by taking the matrix element of this operator between 
the reference state | ti?,)- The full Augmented Space 
Theorem is a generalization of this for functions of many 
independent random variables {tir}. 
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APPENDIX B: OPERATIONS IN AUGMENTED 
RECIPROCAL SPACE 

The main operations of the effective dynamical ma- 
trix on a general state |k £g> {c}) = ||{c}) in the aug- 



mented reciprocal space are given below, where {c} = 
{Ri, i?2i ■ ■ ■ i Rc\ indicates the cardinality sequence with 
| |) configuration at R±, R2, ■ ■ ■ , R c sites. This is required 
to implement the recursion procedure : 



J 



(£ Tl l ®P^J II {C}> 

EE f ® T w) 

. R X ) 



{C}) = 



J2 P l R® p n)^Y, exp(ik-R')\R',{C}) 
{C}) S (R G {C}) (i?o is any reference site) 



{C ± Ro}) 



Z^T^R+x) 4=E exp(*.J?OI^{C}) 
, R x / V R' 



— s(k) || {C — x}) (x is a lattice vector) 
where s(k) = J2 X ex P (— *k • x) 



EE {pr+p 1 r + x)® t R ,R+x 

\ R x 



EE ( 



?r +T£ x )®T R , R+> 



R x 



EE (prPr+x)® t ^+x 

. R x 



EE (T^ +X )®T R , R+X 



(pR + x r R+PR r i 



Tl 

R+x 



R,R+X 



We note that all operations involve only manipulations 
of the configuration part of the state ( i.e. manipula- 
tions of the cardinality sequence only ). The operation of 
the effective dynamical matrix thus entirely takes place 
in the configuration space and the calculation does not 
involve the real space Tt at all. This is an enormous 
simplification over the standard augmented space recur- 
sion described earlier |19|. where the entire reduced real 
space part as well as the configuration part was involved 



{C}) = a (k) || {C - X }> S (Ro e {C - x}) + S(R Q + xe {C - x }) 



|| {C}) = a (k) || {C X } ± Ro) + II {C - x} ± (Ro + x)) 



{C}) = a(k) || {C - x}) [5 (i?o e {C - x}) 5 (Ro + x e {C - X }) 



{C}) = s(k) || {C- X }±Ro±(Ro + x)) 



|| {C}) = a(k) || {C- X }± Ro) S ((Ro + x)e{C- x }) + 



|| {C - x} ± (Ro + x)) S (Ro e {C - x}) 



(Bl) 



r 



in the recursion process. Since one can efficiently store 
the configurations in bits of words so now the calculation 
becomes much simpler. These operations finally involve 
simple bit manipulation routines. 

It is interesting to note that the second operation in 
the above list creates a new configuration. In the next 
step of recursion the third operation translates the en- 
tire operation by lattice translations {x}- The cluster of 
configurations thus 'travel' across the lattice as recursion 
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proceeds. 

APPENDIX C: TERMINATORS 

The recursive calculation described earlier gives rise to 
a set of continued fraction coefficients {a ra ,6„}. In any 
practical calculation we can go only upto a finite number 
of steps, consistent with our computational process. In 
case the coefficients converge, i.e. if \a n — a\ < e, |6„ — b\ < 
e for n > N, we may replace {a n ,b n } by {a, b} for all 
n > N . In that case the asymptotic part of the contiued 
fraction may be analytically summed to obtain : 

T{E) = (1/2) (E-a- V(^-a) 2 - 46 2 ) 
I 



which gives a continuous spectrum a — 2b < E < a + 2b. 
Since the terminator coefficients are related to the band 
edges and widths, a sensible criterion for the choice of 
these asymptotic coefficients is necessary, so as not to give 
arise to spurious structures in our calculations. Beer and 
Pettifor [23] suggest a sensible criterion : given a finite 
number of coefficients, we must choose {a, b} in such a 
way so as to give, for this set of coefficients, the minimum 
bandwidth consistent with no loss of spectral weight from 
the band. This criterion is easily translated into math- 
ematical terms. The delta functions that would carry 
weight out of the band must then be situated exactly 
at the band edges. We thus demand that the continued 
fraction diverge simulatneously at both the top and the 
bottom of the band. 

At the band edges : T(a ± 26) = ±6 so, 



< G(a ± 2b) »= 



b\/A 



±b--(a 2 -a) - 



&§/4 

b%/2 
±b — (ajy — a) 



For a given a, the (N+l) eigenvalues of the finite tridi- 
agonal matrix : 



ai - 
¥2 




\b 2 



|(a 2 - a) 



±63 



V 



V2 



>N 



V2 

(a N 



>N 



a) J 



are values at which the Green function diverges. The 
maximum and minimum of this set of eigenvalues are 
those values of b for which spectral weight has just split 
off from the band. Thus our choice of a is that value 
for which the maximum eigenvalue is the largest and the 
minimum the smallest. Since the continued fraction in- 
volves b 2 then, 



sup b„ 

{a} 



c(a c ) 



inf \b, n 

{a} 



With this choice the terminator T(E) has all the Her- 
glotz properties required. Luchini and Nex 22] further 
modified this by replacing the "butt joining" {a n ,b n } to 
a, b by a smooth linear interpolation : 



,b n 



a n , b n n <n 1 

) (a n (N -n) + a(n-ni)]/(N - m) n x < n < N 
{b n (N -n) + b(n-ni)]/(N -ni) m<n<N 
a,b n > N 

They argued that most of the possible spurious 
structures are removed by such interpolation. In our 
work we have used these two ideas to estimate the 
terminator. 
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